# -*- coding: utf-8 -*-
"""
Created on Mon Nov 15 16:05:10 2021

@author: yhl
"""

"""MyProblem.py"""

import numpy as np
import geatpy as ea
from iapws import IAPWS97



class MVSG (ea.Problem):  # 继承Problem父类
    def __init__(self):
        name = 'MVSG'  # 初始化name（函数名称，可以随意设置）
        M = 1  # 初始化M（目标维数）
        maxormins = [1]  # 初始化maxormins（目标最小最大化标记列表，1：最小化该目标；-1：最大化该目标）
        Dim = 7  # 初始化Dim（决策变量维数）
        varTypes = [0] * Dim  # 初始化varTypes（决策变量的类型，元素为0表示对应的变量是连续的；1表示是离散的）
        lb = [15,  1.2 ,  3,  14,  290.0,  260.0,  3.2] #决策变量下界
        ub = [27,  1.5 ,  5,  17,  320.0,  290.0,  5.2] #决策变量上界
        lbin = [1,1,1,1,1,1,1] #决策变量下边界
        ubin = [1,1,1,1,1,1,1] #决策变量上边界
        # 调用父类构造方法完成实例化
        ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)

    def aimFunc(self, pop):  # 目标函数
        Vars = pop.Phen  # 得到决策变量矩阵
        do = Vars[:, [0]]  # 取出第一列得到所有个体的x1组成的列向量
        xs = Vars[:, [1]]  # 取出第二列得到所有个体的x2组成的列向量
        u1 = Vars[:, [2]]  # 取出第三列得到所有个体的x3组成的列向量
        P1 = Vars[:, [3]]  # 取出第四列得到所有个体的x4组成的列向量
        t1_i = Vars[:, [4]]  # 取出第五列得到所有个体的x5组成的列向量
        t1_ii = Vars[:, [5]]  # 取出第六列得到所有个体的x6组成的列向量
        P2 = Vars[:, [6]]  # 取出第七列得到所有个体的x7组成的列向量
        # 计算目标函数值，赋值给pop种群对象的ObjV属性
        #########################
        # 待优化参数
        #        P1=15 #MPa 一回路定性压力
        #        t1_i=297.705
        #        t1_ii=262.295
        #        u1=5
        #        do=19.05
        #        xs=1.4
        #        P2=3.6 #MPa
        #########################
        # 初始已知热力参数
        Pt = 150000  # kW
        t1 = (t1_i + t1_ii) / 2  # 一回路定性温度
        P2f = 5.2  # MPa
        tf = 105
        ηsg = 0.99
        x = 0.999
        Cs = 0.011
        λw = 18.4
        Rf = 0.0000088
        Cwl = 0.013
        np.pi = 3.141592653
        g = 9.8
        jc = 200
        σ_x = 183
        σ_xf = 146
        σ_s = 143
        bh = P1 * do / (1.5 * σ_x - 0.6 * P1)  # 壁厚
        di = do - 2 * bh
        t = xs * do
        # 热平衡计算
        # 工质状态
        h1_i = np.zeros((do.shape[0], do.shape[1]))
        h1_ii = np.zeros((do.shape[0], do.shape[1]))
        ρ1_i = np.zeros((do.shape[0], do.shape[1]))
        ρ1_ii = np.zeros((do.shape[0], do.shape[1]))

        h1 = np.zeros((do.shape[0], do.shape[1]))
        v1 = np.zeros((do.shape[0], do.shape[1]))
        Cpc = np.zeros((do.shape[0], do.shape[1]))
        μc = np.zeros((do.shape[0], do.shape[1]))
        υc = np.zeros((do.shape[0], do.shape[1]))
        λc = np.zeros((do.shape[0], do.shape[1]))
        Prc = np.zeros((do.shape[0], do.shape[1]))

        ts = np.zeros((do.shape[0], do.shape[1]))
        hs = np.zeros((do.shape[0], do.shape[1]))
        Prs = np.zeros((do.shape[0], do.shape[1]))
        μs = np.zeros((do.shape[0], do.shape[1]))
        ρs = np.zeros((do.shape[0], do.shape[1]))
        ρs_i = np.zeros((do.shape[0], do.shape[1]))
        Cps_i = np.zeros((do.shape[0], do.shape[1]))
        h2_ii = np.zeros((do.shape[0], do.shape[1]))
        hf = np.zeros((do.shape[0], do.shape[1]))
        vf = np.zeros((do.shape[0], do.shape[1]))
        Cpf = np.zeros((do.shape[0], do.shape[1]))

        for i in range(do.shape[0]):
            for j in range(do.shape[1]):
                # 一回路冷却剂进出口参数
                h1_i[i, j] = IAPWS97(P=P1[i, j], T=t1_i[i, j] + 273.15).h
                h1_ii[i, j] = IAPWS97(P=P1[i, j], T=t1_ii[i, j] + 273.15).h
                ρ1_i[i, j] = IAPWS97(P=P1[i, j], T=t1_i[i, j] + 273.15).rho
                ρ1_ii[i, j] = IAPWS97(P=P1[i, j], T=t1_ii[i, j] + 273.15).rho

                # 一回路定性参数下的水热力参数
                h1[i, j] = IAPWS97(P=P1[i, j], T=t1[i, j] + 273.15).h
                v1[i, j] = IAPWS97(P=P1[i, j], T=t1[i, j] + 273.15).v
                Cpc[i, j] = IAPWS97(P=P1[i, j], T=t1[i, j] + 273.15).cp
                μc[i, j] = IAPWS97(P=P1[i, j], T=t1[i, j] + 273.15).mu
                υc[i, j] = IAPWS97(P=P1[i, j], T=t1[i, j] + 273.15).nu
                λc[i, j] = IAPWS97(P=P1[i, j], T=t1[i, j] + 273.15).k
                Prc[i, j] = IAPWS97(P=P1[i, j], T=t1[i, j] + 273.15).Prandt

                # 二回路水热力参数
                ts[i, j] = IAPWS97(P=P2[i, j], x=0).T - 273.15
                hs[i, j] = IAPWS97(P=P2[i, j], x=0).h
                Prs[i, j] = IAPWS97(P=P2[i, j], x=0).Prandt
                μs[i, j] = IAPWS97(P=P2[i, j], x=0).mu
                ρs[i, j] = IAPWS97(P=P2[i, j], x=0).rho
                ρs_i[i, j] = IAPWS97(P=P2[i, j], x=1).rho
                Cps_i[i, j] = IAPWS97(P=P2[i, j], x=0).cp
                h2_ii[i, j] = IAPWS97(P=P2[i, j], x=1).h
                hf[i, j] = IAPWS97(P=P2f, T=tf + 273.15).h
                vf[i, j] = IAPWS97(P=P2f, T=tf + 273.15).v
                Cpf[i, j] = IAPWS97(P=P2f, T=tf + 273.15).cp
        G1 = Pt / (h1_i - h1_ii)
        Q = ηsg * Pt
        r = h2_ii - hs
        D = Q / (r * x + (1 + Cs) * (hs - hf))
        Dbd = D * Cs
        Df = Dbd + D
        # 预热段计算
        # 热力参数
        h1 = np.zeros((do.shape[0], do.shape[1]))
        v1 = np.zeros((do.shape[0], do.shape[1]))
        Cpc = np.zeros((do.shape[0], do.shape[1]))
        μc = np.zeros((do.shape[0], do.shape[1]))
        υc = np.zeros((do.shape[0], do.shape[1]))
        λc = np.zeros((do.shape[0], do.shape[1]))
        Prc = np.zeros((do.shape[0], do.shape[1]))
        for i in range(do.shape[0]):
            for j in range(do.shape[1]):
                h1[i, j] = IAPWS97(P=P1[i, j], T=t1[i, j] + 273.15).h
                v1[i, j] = IAPWS97(P=P1[i, j], T=t1[i, j] + 273.15).v
                Cpc[i, j] = IAPWS97(P=P1[i, j], T=t1[i, j] + 273.15).cp
                μc[i, j] = IAPWS97(P=P1[i, j], T=t1[i, j] + 273.15).mu
                υc[i, j] = IAPWS97(P=P1[i, j], T=t1[i, j] + 273.15).nu
                λc[i, j] = IAPWS97(P=P1[i, j], T=t1[i, j] + 273.15).k
                Prc[i, j] = IAPWS97(P=P1[i, j], T=t1[i, j] + 273.15).Prandt

        # 热预热段
        # 调整预热段两端高度相等

        CR_i = 4
        Dr = (CR_i - 1) * D
        r1 = 0.8
        r2 = 1 - r1
        Dre = Dr / 2 + r1 * Df
        Dle = Dr / 2 + r2 * Df
        t_i = np.ones((do.shape[0], do.shape[1])) * 200
        while True:
            Cp_i = np.zeros((do.shape[0], do.shape[1]))
            for i in range(do.shape[0]):
                for j in range(do.shape[1]):
                    Cp_i[i, j] = IAPWS97(P=P2[i, j], T=t_i[i, j] + 273.15).cp
            _Cp1_i = (Cps_i + Cp_i) / 2
            _Cp2_i = (Cpf + Cp_i) / 2
            _t_i = (Dr * _Cp1_i * ts + Df * _Cp2_i * tf) / (Dr * _Cp1_i + Df * _Cp2_i)
            if np.any((_t_i - t_i) >= 0.01) == 1:
                t_i = t_i - (t_i - _t_i) / 10
            else:
                break

        _h = np.zeros((do.shape[0], do.shape[1]))
        for i in range(do.shape[0]):
            for j in range(do.shape[1]):
                _h[i, j] = IAPWS97(P=P2[i, j], T=t_i[i, j] + 273.15).h
        delta_hy_i = hs - _h
        Q1 = delta_hy_i * (Df + Dr)
        # 预热段冷却剂
        h1b_i = h1_i - (Dre * delta_hy_i / G1)
        t1b_i = np.zeros((do.shape[0], do.shape[1]))
        for i in range(do.shape[0]):
            for j in range(do.shape[1]):
                t1b_i[i, j] = IAPWS97(P=P1[i, j], h=h1b_i[i, j]).T - 273.15
        t1y_i = (t1_i + t1b_i) / 2
        υ1y_i = np.zeros((do.shape[0], do.shape[1]))
        λ1y_i = np.zeros((do.shape[0], do.shape[1]))
        Pr1y_i = np.zeros((do.shape[0], do.shape[1]))
        for i in range(do.shape[0]):
            for j in range(do.shape[1]):
                υ1y_i[i, j] = IAPWS97(P=P1[i, j], T=t1y_i[i, j] + 273.15).nu
                λ1y_i[i, j] = IAPWS97(P=P1[i, j], T=t1y_i[i, j] + 273.15).k
                Pr1y_i[i, j] = IAPWS97(P=P1[i, j], T=t1y_i[i, j] + 273.15).Prandt
        # 预热段二回路水
        ty_i = (t_i + ts) / 2
        Cp2y_i = np.zeros((do.shape[0], do.shape[1]))
        μ2y_i = np.zeros((do.shape[0], do.shape[1]))
        υ2y_i = np.zeros((do.shape[0], do.shape[1]))
        λ2y_i = np.zeros((do.shape[0], do.shape[1]))
        for i in range(do.shape[0]):
            for j in range(do.shape[1]):
                Cp2y_i[i, j] = IAPWS97(P=P2[i, j], T=ty_i[i, j] + 273.15).cp
                μ2y_i[i, j] = IAPWS97(P=P2[i, j], T=ty_i[i, j] + 273.15).mu
                υ2y_i[i, j] = IAPWS97(P=P2[i, j], T=ty_i[i, j] + 273.15).nu
                λ2y_i[i, j] = IAPWS97(P=P2[i, j], T=ty_i[i, j] + 273.15).k
        Re1f_i = u1 * di / 1000 / υ1y_i
        α1y_i = 0.023 * Re1f_i ** 0.8 * Pr1y_i ** 0.4 * λ1y_i / di * 1000
        Rw = (do / 1000) / (2 * λw) * np.log(do / di)

        # 假定校核热预热段热流密度
        q1 = 31425
        while True:
            twy_i = t1y_i - q1 * (1 / α1y_i * do / di + Rw + Rf)
            μwy_i = np.zeros((do.shape[0], do.shape[1]))
            for i in range(do.shape[0]):
                for j in range(do.shape[1]):
                    μwy_i[i, j] = IAPWS97(P=P2[i, j], T=twy_i[i, j] + 273.15).mu
            de = (4 * t ** 2 - np.pi * do ** 2) / np.pi / do  # 单位：mm
            uf = 0.5  # 选定
            Re2y_i = uf * de / 1000 / υ2y_i
            α2y = jc * λ2y_i / de * 1000 * (1000 * Cp2y_i * μ2y_i / λ2y_i) ** (1 / 3) * (μ2y_i / μwy_i) ** 0.14
            ky_i = 1 / (1 / α1y_i * do / di + Rw + Rf + 1 / α2y)
            delta_tmax_y = t1_i - t_i
            delta_tmin_y = t1b_i - ts
            delta_tm_y = (delta_tmax_y - delta_tmin_y) / np.log(delta_tmax_y / delta_tmin_y)
            q11 = ky_i * delta_tm_y
            if np.any(abs(q1 - q11) > 10) == 1:
                q1 = q11
            else:
                break
        #        print("q1",q1)
        # 冷预热段
        # 预热段冷却剂
        h1b_ii = h1_ii + Dle * delta_hy_i / G1
        t1b_ii = np.zeros((do.shape[0], do.shape[1]))
        for i in range(do.shape[0]):
            for j in range(do.shape[1]):
                t1b_ii[i, j] = IAPWS97(P=P1[i, j], h=h1b_ii[i, j]).T - 273.15
        t1y_ii = (t1_ii + t1b_ii) / 2
        λ1y_ii = np.zeros((do.shape[0], do.shape[1]))
        υ1y_ii = np.zeros((do.shape[0], do.shape[1]))
        Pr1y_ii = np.zeros((do.shape[0], do.shape[1]))
        for i in range(do.shape[0]):
            for j in range(do.shape[1]):
                λ1y_ii[i, j] = IAPWS97(P=P1[i, j], T=t1y_ii[i, j] + 273.15).k
                Pr1y_ii[i, j] = IAPWS97(P=P1[i, j], T=t1y_ii[i, j] + 273.15).Prandt
                υ1y_ii[i, j] = IAPWS97(P=P1[i, j], T=t1y_ii[i, j] + 273.15).nu
        Re1f_ii = u1 * di / 1000 / υ1y_ii
        α1y_ii = 0.023 * Re1f_ii ** 0.8 * Pr1y_ii ** 0.4 * λ1y_ii / di * 1000
        # 假定迭代q2
        q2 = 40000
        for a in range(100):
            twy_ii = t1y_ii - q2 * (1 / α1y_ii * do / di + Rw + Rf)
            μwy_ii = np.zeros((do.shape[0], do.shape[1]))
            for i in range(do.shape[0]):
                for j in range(do.shape[1]):
                    μwy_ii[i, j] = IAPWS97(P=P2[i, j], T=twy_ii[i, j] + 273.15).mu
            α2y_ii = jc * λ2y_i / de * 1000 * (1000 * Cp2y_i * μ2y_i / λ2y_i) ** (1 / 3) * (μ2y_i / μwy_ii) ** 0.14
            ky_ii = 1 / (1 / α1y_ii * do / di + Rw + Rf + 1 / α2y_ii)
            delta_tmax_yy = t1_ii - t_i
            delta_tmin_yy = t1b_ii - ts
            delta_tm_yy = (delta_tmax_yy - delta_tmin_yy) / np.log(delta_tmax_yy / delta_tmin_yy)
            q21 = ky_ii * delta_tm_yy
            if np.any(abs(q2 - q21) > 10) == 1:
                q2 = q2 - (q2 - q21) / 10
            else:
                break
        #        print("q2",q2)
        # 预热段总面积
        Cy = 1.08
        Fy = Cy * 2 * Q1 * 1000 / (q1 + q2)
        _Fy = Fy / 2
        # 预热段换热量
        Qy_i = q1 * _Fy / 1000
        Qy_ii = q2 * _Fy / 1000
        # 沸腾段
        Q2 = Q - Q1
        t1b = (t1b_i + t1b_ii) / 2
        υ1b = np.zeros((do.shape[0], do.shape[1]))
        λ1b = np.zeros((do.shape[0], do.shape[1]))
        Pr1b = np.zeros((do.shape[0], do.shape[1]))
        for i in range(do.shape[0]):
            for j in range(do.shape[1]):
                υ1b[i, j] = IAPWS97(P=P1[i, j], T=t1b[i, j] + 273.15).nu
                λ1b[i, j] = IAPWS97(P=P1[i, j], T=t1b[i, j] + 273.15).k
                Pr1b[i, j] = IAPWS97(P=P1[i, j], T=t1b[i, j] + 273.15).Prandt
        Re1b = u1 * di / 1000 / υ1b
        α1b = 0.023 * Re1b ** 0.8 * Pr1b ** 0.3 * λ1b / di * 1000

        delta_tmax_b = t1b_i - ts
        delta_tmin_b = t1b_ii - ts
        delta_tm_b = (delta_tmax_b - delta_tmin_b) / np.log(delta_tmax_b / delta_tmin_b)
        # 假定迭代热流密度qb
        qb = 100000
        while True:
            σ = 0.027410432  # 没查表
            α2b = (Cps_i / r / Cwl / Prs) * (μs * r * 1000 * (g * (ρs - ρs_i) / σ) ** 0.5) ** 0.33 * (qb ** 0.67)
            kb = 1 / (1 / α1b * do / di + Rw + Rf + 1 / α2b)
            qb1 = kb * delta_tm_b
            if np.any(abs(qb - qb1) > 10) == 1:
                qb = qb - (qb - qb1) / 10
            else:
                break
        #        print("qb",qb)
        Fb = Q2 * 1000 / qb
        Cb = 1.08
        F = Cb * Fb + Fy

        # 结构计算
        # U型管设计
        R = 2 * t
        A = G1 * v1 / u1
        a = 1 / 4 * np.pi * di ** 2
        n_ii = A / a * 1000000
        C = 1.08
        n = np.ceil(C * n_ii)

        # 管板计算
        nk = 2 * n
        Dt = 10000
        while True:
            n_i = (Dt - do) / t + 1
            tn = 2 * R
            Ad = n_i * t * (tn - t)
            At = 2 * n * t ** 2 + Ad
            Dt_i = (4 * At / np.pi) ** 0.5
            if np.any((Dt - Dt_i) > 1) == 1:
                Dt = Dt - (Dt - Dt_i) / 10
            else:
                break
        Dt = np.ceil(Dt)
        e = 30
        Dgb = Dt + e
        σgb = 131.3
        Sdh_i = 5
        Sgb = np.ceil(0.5 * 1.04 * Dgb * (P1 / σgb) ** 0.5 + Sdh_i)
        #        print("n",n)
        # 管束计算
        Lz = np.ceil(F / np.pi / do * 1000)
        Dj = 4 * t
        Rmax = Dt / 2
        R = (Rmax + 0.5 * Dj) / 2
        Lw = n * np.pi * R / 1000
        Lzh = Lz - Lw
        Hw = Rmax / 1000
        Hz = Lzh / 2 / n
        Htb = Hw + Hz
        Dtb = 2 * Hw
        Cgj = Htb / Dtb
        l = Lz / n + 2 * Sgb / 1000
        lmax = np.pi * Rmax / 1000 + 2 * Hz + 2 * Sgb / 1000
        lmin = np.pi * Dj / 2 / 1000 + 2 * Hz + 2 * Sgb / 1000

        # 管束套筒、下封头及上下筒体
        P1d = P1
        δt = 0.02
        Dwi = (Dtb + 2 * δt) * 1000
        δ = 12
        Dwo = Dwi + 2 * δ
        B = 88
        Dix = Dwo + 2 * B
        # Sx=84
        ############################
        Sx = P2 * Dix / 2 / (2 * σ_x - 0.6 * P2)
        Sx = np.ceil(Sx)
        ############################
        Dox = Dix + 2 * Sx
        Sl = Sx + 10
        Do_xf = Dox + 2 * Sl
        ##############################
        Sxf = P1d * Do_xf / 2 / (2 * σ_xf - 0.6 * P1d)
        Sxf = np.ceil(Sxf)
        ##############################
        Di_xf = Do_xf - 2 * Sxf
        # 各主要部分内经参数
        d1i = 374
        u1_ii = 9.94

        # 水力计算
        # 一回路流动阻力
        u1_i = 1.05 * u1
        Re = u1_i * di / 1000 / υc
        f = 0.184 * Re ** (-0.2)
        delta_Pf = f * l / di * 1000 * u1_i ** 2 / 0.989 / v1  # 单位：Pa
        # print("delta_Pf",delta_Pf)
        Fc = np.pi * (Di_xf / 1000) ** 2 / 8
        A1 = np.pi * (d1i / 1000) ** 2 / 4
        ξ1 = (1 - A1 / Fc) ** 2
        delta_P1 = ξ1 * ρ1_i * u1_ii ** 2 / 2
        ξ2 = 0.9
        delta_P2 = ξ2 * ρ1_i * u1_ii ** 2 / 2
        A_i = A / 1.05
        ξ3 = (1 - A_i / Fc) ** 2
        delta_P3 = ξ3 * ρ1_i * u1_i ** 2 / 2
        ξ4 = 0.5
        delta_P4 = ξ4 * ρ1_i * u1_i ** 2 / 2
        ξ5 = ξ3
        delta_P5 = ξ5 * ρ1_ii * u1_i ** 2 / 2
        ξ6 = 0.9
        delta_P6 = ξ6 * ρ1_ii * u1_ii ** 2 / 2
        ξ7 = ξ1
        delta_P7 = ξ7 * ρ1_ii * u1_ii ** 2 / 2
        delta_Pl = delta_P1 + delta_P2 + delta_P3 + delta_P4 + delta_P5 + delta_P6 + delta_P7
        delta_P = 1.1 * (delta_Pf + delta_Pl)
        # 二回路下降空间阻力
        vd = np.zeros((do.shape[0], do.shape[1]))
        for i in range(do.shape[0]):
            for j in range(do.shape[1]):
                vd[i, j] = IAPWS97(P=P1[i, j], T=t_i[i, j] + 273.15).v
        ρd = 1 / vd
        ξin = 0.95
        ξout = 0.95
        ξt = 0.95
        ξf = 0.95
        Hd = Htb + 1.5  # 下降空间高度 设定
        De = Dix - Dwo
        delta = 0.15
        fd = (1.74 + 2 * np.log10(De / 2 / delta)) ** (-2)
        vs = np.zeros((do.shape[0], do.shape[1]))
        v2_ii = np.zeros((do.shape[0], do.shape[1]))
        υl = np.zeros((do.shape[0], do.shape[1]))
        υs = np.zeros((do.shape[0], do.shape[1]))
        for i in range(do.shape[0]):
            for j in range(do.shape[1]):
                vs[i, j] = IAPWS97(P=P2[i, j], x=0).v
        ρs = 1 / vs
        for i in range(do.shape[0]):
            for j in range(do.shape[1]):
                v2_ii[i, j] = IAPWS97(P=P2[i, j], x=1).v
        ρ2_ii = 1 / v2_ii
        Fu = np.pi / 4 * ((Dwi / 1000) ** 2 - 2 * n * (do / 1000) ** 2)
        de_2 = 4 * Fu * 1000 / np.pi / (Dwi + 2 * n * do)  # 原为de 单位：m
        for i in range(do.shape[0]):
            for j in range(do.shape[1]):
                υl[i, j] = IAPWS97(P=P2[i, j], x=0).nu
                υs[i, j] = IAPWS97(P=P2[i, j], x=1).nu

        CR = np.ones((do.shape[0], do.shape[1])) * 2.9
        for a in range(220):
            #        while True:
            # 下降空间阻力
            Fd = np.pi / 4 * ((Dix / 1000) ** 2 - (Dwo / 1000) ** 2)
            ud = CR * D * vd / Fd
            delta_Pd = (fd * Hd / De * 1000 + ξin + ξout + ξt + ξf) * ρd * ud ** 2 / 2
            #            print("delta_Pd",delta_Pd)

            # 上升空间阻力
            # 平均摩擦阻力
            uo = CR * D * vs / Fu
            uo2_i = (CR - 1) * D * vs / Fu
            _uo_i = 0.5 * (uo2_i + uo)
            uo_ii = D * v2_ii / Fu
            _uo_ii = 0.5 * uo_ii
            Relo = _uo_i * de_2 / υl
            Reso = _uo_ii * de_2 / υs
            flo = 0.184 * Relo ** (-0.2)
            fso = 0.184 * Reso ** (-0.2)
            delta_Pflo = flo * Hz / de_2 * ρs * _uo_i ** 2 / 2
            delta_Pfso = 1 / 3 * fso * Hz / de_2 * ρ2_ii * _uo_ii ** 2 / 2
            X1 = (delta_Pflo / delta_Pfso) ** 0.5
            C1 = 20
            Φl_21 = 1 + C1 / X1 + 1 / X1 ** 2
            Φs_21 = 1 + C1 * X1 + X1 ** 2
            delta_Pfl = Φl_21 * delta_Pflo
            delta_Pfs = Φs_21 * delta_Pfso
            delta_Pff = 0.5 * (delta_Pfl + delta_Pfs)  # 原delta_Pf
            # print(delta_Pff)

            # 平均局部阻力
            N = 9
            ξ11 = 3
            delta_Pllo1 = N * ξ11 * ρs * _uo_i ** 2 / 2
            delta_Plso1 = 1 / 3 * N * ξ11 * ρ2_ii * _uo_ii ** 2 / 2
            X2 = (delta_Pllo1 / delta_Plso1) ** 0.5
            Pc = 22.06
            ZR = 1 / (0.19 + 0.92 * P2 / Pc)
            K = ZR + 1 / ZR
            Φl_22 = 1 + K / X2 + 1 / X2 ** 2
            Φs_22 = 1 + K * X2 + X2 ** 2
            delta_Pll = Φl_22 * delta_Pllo1
            delta_Pls = Φs_22 * delta_Plso1
            delta_Pl = 0.5 * (delta_Pll + delta_Pls)
            # print(delta_Pl)

            # 弯头区局部阻力
            Djmax = 2 * Hw - do / 1000 / 2
            ys = 0.2122 * Djmax
            N_i = ys / t * 1000 - 1
            # print(N_i)
            n_xs = 0.43 + 1.13 / xs
            Relo = uo2_i * de_2 / υl
            Reso = uo_ii * de_2 / υs
            flo = 4 * (0.044 + 0.08 * xs / (xs - 1) ** n_xs) * Relo ** (-0.15)
            fso = 4 * (0.044 + 0.08 * xs / (xs - 1) ** n_xs) * Reso ** (-0.15)  # 乘4后与算例相同
            delta_Pblo = N_i * flo * ρs * uo2_i ** 2 / 2
            delta_Pbso = 1 / 3 * N_i * fso * ρ2_ii * uo_ii ** 2 / 2
            X3 = (delta_Pblo / delta_Pbso) ** 0.5
            C3 = 20
            Φl_23 = 1 + C3 / X3 + 1 / X3 ** 2
            Φs_23 = 1 + C3 * X3 + X3 ** 2
            delta_Pbl = Φl_23 * delta_Pblo
            delta_Pbs = Φs_23 * delta_Pbso
            delta_Pb = 0.5 * (delta_Pbl + delta_Pbs)

            # 加速度阻力
            x2 = 1 / CR
            β2 = x2 / ρ2_ii / (x2 / ρ2_ii + (1 - x2) / ρs)
            C = 0.833 + 0.05 * np.log(P2)
            φ2 = C * β2
            G_ = uo * ρs
            delta_Pa = G_ ** 2 * ((1 - x2) ** 2 / ρs / (1 - φ2) + x2 ** 2 / ρ2_ii / φ2 - 1 / ρs)
            # print(delta_Pa)

            # 流量分配孔板局部阻力
            ξh = 3
            delta_Ph = ξh * ρs * uo ** 2 / 2

            # 上升空间总阻力
            delta_Pr = delta_Pff + delta_Pl + delta_Pb + delta_Pa + delta_Ph

            # 汽水分离器总阻力
            # delta_Ps=2245*CR+5883
            delta_Ps = 2245 * CR + 25880

            # 二回路循环总阻力
            delta_Pz = delta_Pd + delta_Pr + delta_Ps
            #            print("delta_Pz",delta_Pz)
            # 二回路循环压头
            Plow = P2 * 1000000 + ρd * g * Hd  # 单位:Pa
            hsl = np.zeros((do.shape[0], do.shape[1]))
            for i in range(do.shape[0]):
                for j in range(do.shape[1]):
                    hsl[i, j] = IAPWS97(P=Plow[i, j] / 1000000, x=0).h
            delta_h = (hs - hf) * 1000 / CR
            delta_P1 = Plow - P2 * 1000000
            scale = delta_h / delta_P1
            qy_i = q1
            # Dcir=CR*D
            Dre = 0.5 * (CR - 1) * D + r1 * Df
            Hp = ((hs - hf) * 1000 / CR + scale * (ρd * g * Hd - delta_Pd)) / (
                        2 * np.pi * do / 1000 * n * qy_i / Dre + scale * ρd * g)
            Hr1 = Htb - Hp
            Hr2 = Hd - Hr1
            _x1 = x2 / 2
            _β1 = _x1 / ρ2_ii / (_x1 / ρ2_ii + (1 - _x1) / ρs)
            _φ1 = C * _β1
            _φ2 = φ2
            Pm1 = (ρs - ρ2_ii) * g * _φ1 * Hr1
            Pm2 = (ρs - ρ2_ii) * g * _φ2 * Hr2
            Pm = Pm1 + Pm2
            #            print("Pm",Pm)
            for i in range(do.shape[0]):
                for j in range(do.shape[1]):
                    if np.any((Pm[i, j] - delta_Pz[i, j]) > 100) == 1:
                        CR[i, j] = CR[i, j] + 0.01
                    else:
                        break
        # print("CR",CR)
        # 上筒体参数
        Dis = Dix + 860  # 设定
        Ss = P2 * Dis / 2 / (2 * σ_s - 0.6 * P2)
        Ss = np.ceil(Ss)
        Dos = Dis + 2 * Ss
        HH = 2.1  # m 上筒体高度 假设
        # 锥形筒体高度与半顶角度
        α = 60
        S = 150
        E = 1
        Dc = (Dix + Dis) / 2
        t3 = P2 * Dc / 1000 / (2 * np.cos(α / 90 * np.pi / 2) * (S * E - 0.6 * P2))
        # 上封头参数
        Di = Dis / 1000
        t2 = P2 * Di / (2 * S - 0.2 * P2)
        Dii = Dos + t2 * 2
        h1 = Dii / 4
        Hs = 1 / 3 * Htb  # m 上封头总高度 假设
        h = Hs - h1 / 1000
        h1_i = (Dii / 1000 + 2 * t2) / 4
        # 重量计算
        ρ1 = 7.85  # Mg/m3 管板、封头：碳钢
        ρ2 = 8.19  # Mg/m3 管束：因科镍690
        # 下封头重量
        W1 = 2 / 3 * np.pi * ρ1 * ((Do_xf / 2000) ** 3 - (Di_xf / 2000) ** 3)
        W1 = np.ceil(W1 * 1000)
        # 管板重量
        W2 = np.pi / 4 * ((Dgb / 1000) ** 2 - nk * (do / 1000) ** 2) * ρ1 * Sgb / 1000
        W2 = np.ceil(W2 * 1000)
        # 管束重量
        W3 = np.pi / 4 * ((do / 1000) ** 2 - (di / 1000) ** 2) * Lz * ρ2
        W3 = np.ceil(W3 * 1000)
        # 下筒体重量
        W4 = np.pi / 4 * ((Dox / 1000) ** 2 - (Dix / 1000) ** 2) * ρ1 * Htb
        W4 = np.ceil(W4 * 1000)
        # 锥形筒体重量
        W5 = np.pi / 24 * ((Dis / 1000 + 2 * t3) ** 3 - (Dix / 1000 + 2 * t3) ** 3 - (Dis / 1000) ** 3 + (
                    Dix / 1000) ** 3) * ρ1 / np.tan(α / 90 * np.pi / 2)
        W5 = np.ceil(W5 * 1000)
        # 上筒体重量
        W6 = np.pi / 4 * ((Dos / 1000) ** 2 - (Dis / 1000) ** 2) * HH * ρ1
        W6 = np.ceil(W6 * 1000)
        # 上封头重量
        W7 = (np.pi / 24 * ((Dii / 1000 + 2 * t2) ** 3 - (Dii / 1000) ** 3) + np.pi / 4 * (
                    (Dii / 1000 + 2 * t2) ** 2 - (Dii / 1000) ** 2) * h) * ρ1
        W7 = np.ceil(W7 * 1000)
        # 汽水分离器重量
        W8 = 4000
        # 总重量
        W_total = (W1 + W2 + W3 + W4 + W5 + W6 + W7 + W8) / 1000
        # print("W1",W1,"W2",W2,"W3",W3,"W4",W4,"W5",W5,"W6",W6,"W7",W7,"W8",W8)
        # 总体积
        V1 = 0.5 * 4 / 3 * np.pi * (Do_xf / 1000 / 2) ** 3
        V4 = np.pi / 4 * (Dox / 1000) ** 2 * Htb
        V5 = np.pi / 24 * ((Dos / 1000) ** 3 - (Dox / 1000) ** 3) / np.tan(α / 90 * np.pi / 2)
        V6 = np.pi / 4 * (Dos / 1000) ** 2 * HH
        V7 = np.pi / 4 * (Di + 2 * t2) ** 3 + np.pi / 4 * (Di + 2 * t2) ** 2 * h
        # print(V1,V4,V5,V6,V7)
        V_total = V1 + V4 + V5 + V6 + V7
        a1 = 0.5
        a2 = 1 - a1
        H_total = a1 * W_total / 22.52 + a2 * V_total / 47.40

        #        采用可行性法则处理约束，生成种群个体违反约束程度矩阵
        for i in range(pop.ObjV.shape[0]):
            for j in range(pop.ObjV.shape[1]):
                if ((Cgj[i, j] < 2 or Cgj[i, j] > 5) or Hr2[i, j] < 0 or D[i, j] < 55 or n[i, j] > 1500 \
                        or lmax[i, j] > 15 or lmin[i, j] < 7 or Htb[i, j] > 8 or Dt[i, j] > 2000 \
                        or delta_P[i, j] > 500000 or (CR[i, j] < 3 or CR[i, j] > 5) or ud[i, j] < 0.3):
                    print("i,j:", i, j)
                    print("pop.ObjV[i,j]:", pop.ObjV[i, j])
                    pop.ObjV[i, j] = max(pop.ObjV)
                    #                    pop.ObjV[i,j]=1000
                    print("max(pop.ObjV):", max(pop.ObjV))
        print("pop.ObjV:", (pop.ObjV))
        return pop.ObjV


        ##############################################################################################################
            # print("P1",P1,"t1_i",t1_i,"t1_ii",t1_ii,"u1",u1,"do",do,"xs",xs,"P2",P2)
            # print("重量："+str(W_total),"体积："+str(V_total),"循环倍率："+str(CR),"总传热面积："+str(F),"传热管根数："+str(n),"r1："+str(r1))
            # print("重量："+str(W_total/232.418),"体积："+str(V_total/103.333),"循环倍率："+str(CR/4.04),"总传热面积："+str(F/6602.5),"传热管根数："+str(n/5438),\
            #    "一回路压力："+str(P1/15.5),"二回路压力："+str(P2/6.89),"进口温度："+str(t1_i/329.8),"出口温度："+str(t1_ii/292.4),\
            #    "传热管外径："+str(do/19.05),"节径比："+str(xs/1.4),"冷却剂流速："+str(u1/5))
            # print("重量循环结束")
            # 采用可行性法则处理约束，生成种群个体违反约束程度矩阵


    def calReferObjV(self):  # 设定目标数参考值（本问题目标函数参考值设定为理论最优值）
        referenceObjV = np.array([[2.5]])
        return referenceObjV



